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ABSTRACT 

Digital  filtering  is  the  process  of  spectrum  shaping  using  digital  components 
as  the  basic  elements.  Increasing  speed  and  decreasing  size  and  cost  of  digital  com¬ 
ponents  make  it  likely  that  digital  filtering,  already  used  extensively  in  the  computer 
simulation  of  analog  filters,  will  perform,  in  real-time  devices,  the  functions  which 
are  now  performed  almost  exclusively  by  analog  components.  In  this  paper,  using 
the  z -transform  calculus,  several  digital  filter  design  techniques  are  reviewed,  and 
new  ones  are  presented.  One  technique  can  be  used  to  design  a  digital  filter  whose 
impulse  response  is  like  that  of  a  given  analog  filter;  another  technique  is  suitable  for 
the  design  of  a  digital  filter  meeting  frequency  response  criteria.  A  third  technique 
yields  digital  filters  with  linear  phase,  specified  frequency  response,  and  controlled 
impulse  response  duration.  The  effect  of  digital  arithmetic  on  the  behavior  of  digital 
filters  is  also  considered. 
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DIGITAL  FILTER  DESIGN  TECHNIQUES 


SECTION  1 
a.  Introduction 

Digital  filtering  is  the  process  of  spectrum  shaping  using  digital  hardware 
as  the  basic  building  block.  Thus  the  aims  of  digital  filtering  are  the  same  as 
those  of  continuous  filtering,  but  the  physical  realization  is  different.  Continuous 
filter  theory  is  based  on  the  mathematics  of  linear  differential  equations;  digital 
filter  theory  is  based  on  the  mathematics  of  linear  difference  equations. 

Interest  in  digital  filtering  derives  from  recent  work  on  computer  simu¬ 
lation  of  speech  communications  systems*’  4’  and  from  problems  arising 
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in  the  processing  of  geophysical  data  .  In  addition,  however,  increasing  speeds 
and  decreasing  costs  of  microelectronic  digital  circuitry  make  possible  the 
conception  of  real  time  digitized  systems  which  perform  the  filtering  operations 
presently  performed  exclusively  by  analog  hardware. 

The  mathematical  basis  of  digital  filter  theory  is  well  known.  In  addition 
to  the  mathematics  literature  on  difference  equations,  a  large  body  of  work  exists 
on  the  subject  of  sampled  data  systems  and  this  work  is  of  direct  use  in  the  design 
of  digital  filters.  Presently,  no  thorough  treatment  exists  of  digital  filter  design 
techniques  and  this  paper  is  a  start  in  that  direction.  Emphasis  will  be  placed  on 
the  frequency  selectivity  of  filters  rather  than  on  the  questions  of  stability  and  time 
response  which  are  usually  of  interest  in  control  systems.  However,  it  should  be 
stressed  that  digital  filters  are  mathematically  equivalent  to  continuous  systems 
with  sampled  data  inputs  and  outputs.  In  the  continuous  case,  the  sampling  interval 
T  usually  depends  on  the  nature  of  the  discrete  input,  such  as  a  pulsed  radar  return. 
The  constraint  on  the  sampling  interval  of  a  digital  filter  could  also  be  caused  by 
the  computation  time  needed  to  execute  the  difference  equations. 

The  basic  mathematical  tool  of  digital  filters  is  the  z  transform  calculus, 

g 

used  by  Hurewicz  to  develop  his  theory  of  pulsed  filters.  We  will  first  briefly 
review  the  z  transform  method  and  then  present  a  number  of  design  techniques  with 
illustrative  examples. 

Real  time  digital  filters  should  have  several  advantages  over  continuous 
filters.  They  should  adhere  to  the  design  so  that  no  tuning  would  be  required.  A 
greater  variety  of  filters  should  be  feasible,  since  no  realizability  problems  (akin 
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to  negative  elements)  arise.  No  special  components  are  needed  to  realize  filters 
with  time  varying  coefficients.  Aggregates  of  digital  filters  should  be  especially 
economic  in  the  very  low  frequency  band  (.  01  cycles  per  second  to  1  cps),  where 
the  size  of  analog  components  becomes  appreciable. 

b.  Difference  Equations  of  Digital  Filters 

An  m  order  linear  difference  equation  may  be  written  as: 

r  m 

y(nT)=£  L  x(nT-iT)-  £  K.  y(nT-  iT)  (1) 

i=o  i=l  1 


The  form  of  (1)  emphasizes  the  iterative  nature  of  the  difference  equation;  given 
the  m  previous  values  of  the  output  y  and  the  r  +  1  most  recent  values  of  the  input 
x,  the  new  output  may  be  computed  from  (1).  Physically,  the  input  numbers  are 
samples  of  a  continuous  waveform  and  real  time  digital  filtering  consists  of  per¬ 
forming  the  iteration  (1)  for  each  arrival  of  a  new  input  sample.  Design  of  the 
filter  consists  of  finding  the  constants  IC  and  to  fulfill  a  given  filtering  require¬ 
ment.  Real  time  filtering  implies  that  the  execution  time  of  the  computer  program 
for  computing  the  right  side  of  (1)  be  less  than  T,  the  sampling  interval.  Simula¬ 
tion  programs  in  which  the  computations  are  not  carried  out  in  real  time  are 
understood  to  be  simulations  of  the  real  time  digital  filters  and  the  original  sampling 
interval  T  of  the  continuous  input  data  remains  the  reference  parameter.  It  should 
be  realized  that  quantization  as  well  as  sampling  is  performed  on  the  continuous 
data  before  entry  into  the  program  and  also,  that  the  multiplications  indicated  in 
(1)  contain  an  inherent  round  off  error  caused  by  finite  register  lengths,  but  for  the 
most  part  these  effects  will  be  neglected.  The  effect  of  quantization  is  an  important 
consideration  and  we  will  study  it,  to  some  extent,  theoretically. 

A  useful  pictorial  representation  of  (1)  consisting  of  unit  delays  (of  time  T), 
adders  and  multipliers  is  shown  in  Fig.  1. 

Many  practical  designs  utilize  arrangements  of  simple  (1st  and  2nd  order) 
difference  equations.  Arrangements  can  be  serial,  parallel  or  combinations  of 
both.  For  example,  the  two  equations: 


y  1(nT)  =  Kuyi(nT  -  T)  +  K^y^nT  -  2T)  +  LR  x(nT) 


y2(nT)  =  K21y2(nT  -  T)  +  K22y2(nT  -2T)  +  L^y^nT) 


(2) 


2 


constitute  a  serial  arrangement  whereby  the  output  y  ^  of  the  first  equation  is 
used  as  input  to  the  second  equation,  as  seen  in  Fig.  2.  Similarly  the  equations: 


yi(nT)  =  K11y1(nT-T)  +  K12y1(nT-2T)+Lnx(nT) 
y2(nT)  =  K21y2(nT  -  T)  +  K22y2(nT  -  2T)  +  L21  x(nT)  S 
y3(nT)  =  y1(nT)  +  y2(nT) 


(3) 


shown  in  Fig.  3,  constitute  a  parallel  arrangement.  The  sets  of  equations  (2) 

th. 

and  (3)  can  be  rewritten  as  single  4Ln  order  difference  equations;  however,  the 
significant  properties  of  the  digital  filters  are  more  readily  understood  via 
structures  like  Figs.  2  and  3. 
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Substantial  literature  exists  on  the  problems  of  pre -sampling  filtering 
of  the  input  and  ultimate  filtering  of  the  output  sequence  to  convert  it  to  a  continuous 
wave.  We  shall  restrict  ourselves  to  studying  only  the  sequences. 

c.  The  Z  Transform 

In  studying  physical  systems  it  is  convenient  to  assume  that  the  input  begins 
at  zero  time  and  has  magnitude  zero  beforehand.  Thus  x(nT)  =  0  for  n  <  0  and 
therefore,  in  an  inertial  system  such  as  that  of  equations  (1),  (2)  or  (3),  y(nT)  =  0 
for  n  <  0. 

The  z  transform  of  a  sequence  x(nT)  is  defined  as: 

X(z)  =  f  x(nT)z'n  (4) 

n=o 


For  many  sequences,  the  infinite  sum  of  (4)  can  be  expressed  in  closed  form. 

For  example,  the  z  transform  of  the  sequence  x(nT)  =  0  for  n  <  0,  x(nT)  =  1 

for  nsO  is  X(z)  =  L  z  n  =  - - — T. 

n=o  i  —  z" 1 

The  transformation  variable  z  is,  in  general,  a  complex  variable  and  X(z) 
is  a  function  of  a  complex  variable. 
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The  z  transform  of  a  convergent  sequence  uniquely  defines  that  sequence. 

By  multiplying  both  sides  of  equation  (4)  by  z11  and  integrating  around  any  closed 
curve  in  the  z  plane  which  encloses  all  the  singular  points  (poles)  of  X(z)  and  the 
origin,  we  obtain: 

x(nT)  =  2¥j  f  X(z)  z11" 1  dz  (5) 

where  j  =  /—  1 

The  inverse  z  transform  (5)  explicitly  determines  the  sequence  x(nT)  asso¬ 
ciated  with  a  given  X(z).  For  sequences  which  converge  to  zero  the  unit  circle 
can  be  the  closed  curve  of  integration. 

d.  Solution  of  Difference  Equations  by  z  Transform  Techniques 

th 

The  z  transform  will  now  be  used  to  obtain  an  explicit  solution  of  the  m 
order  linear  difference  equation  (1).  First,  rewrite  (1)  as: 


m  r 

T  K  y(nT-iT)  =  T  L.  x(nT-iT)  (6) 

i=o  i  i=o  1 

where  K  =  1 
o 

Next,  take  the  z  transform  of  (6),  to  obtain: 


m  ®  -n  r 

r  K.  £  y(nT  -  iT)  z  =  T 

i=o  i  n=o  i=o 


L.  T,  x(nT-iT)z"n  (7) 

1  n=o 


Recognizing  that  the  z  transform  of  a  sequence  delayed  by  i  samples  is 
equal  to  the  z  transform  of  the  original  sequence  multiplied  by  z  1,  equation 
(7)  reduces  to: 


m 


-l 


or 


Y(z)  L  K.  z  =  X(z)  D  L.  z 

i=o  1  i=o  1 


•I  ^  z” 1 

Y(z)  =  X(z)  TfT - -  rr  =  X(z)  H(z) 


X  K.  z 
1=0  1 


-1 


(8) 
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Thus  Y(z)  is  explicitly  determined  as  the  product  of  the  input  z  transform 
and  a  system  function  H(z)  which  is  a  rational  fraction  in  z  and  is  a  function  of 
the  constant  coefficients  in  the  original  difference  equation.  By  the  inverse  z 
transform  (5),  the  output  y(nT)  may  be  written  as  the  inverse  transform  of  Y(z). 

If  we  choose  as  the  input  sequence  x(nT)  =  for  ,  then  in  Eq.  (8) 
X(z)  =  1.  Clearly  the  response  to  this  input  is  the  inverse  z  transform  of  H(z). 

Thus  the  sequence  1,  0,  0, .  plays  the  part  in  digital  filter  theory  which  the 

unit  impulse  plays  in  continuous  filter  theory.  We  shall  refer  to  the  inverse  z 
transform  of  H(z)  as  the  impulse  response  of  a  digital  filter. 

We  can  use  Eq.  (8)  to  derive  a  general  representation  of  a  digital  filter 
which  differs  from  that  of  Fig.  1.  We  can  define  an  intermediate  variable  W(z) 
corresponding  to  an  intermediate  sequence  w(nT),  such  that 


W(z)  =  X(z)  * 

m  : 

.£  K.  z 
1=0  1 

Y(z)  =  W(z)  L  L.  z"1 
i=o  i 


(8a) 


The  equations  of  (8a)  result  in  simultaneous  difference  equations 


m 

w(nT)  =  x(nT)  —  .2)^  K.  w(nT  —  iT) 
r 

y(nT)  =  T  L.  w(nT  -  iT)  (8b) 

i=o  1 

which  lead  to  the  structure  shown  in  Fig.  4.  This  structure  requires  fewer  delays 
(less  storage)  and  may,  therefore,  be  preferred  for  some  realizations. 

The  primary  importance  of  the  system  function  H(z)  for  our  purposes  is  in 
its  interpretation  as  a  frequency  selective  function.  Let  us  assume  that  the  input 
is  a  sampled  complex  exponential  wave. 

x(nT)  =  ejnwT  (9) 
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The  solution  to  (6)  for  the  input  (9)  is  also  a  complex  exponential  wave 
which  can  be  represented  as: 


y(nT)  =  F(eJ  )  e 


jncoT 


(10) 


and,  by  substitution  of  (9)  and  (10)  into  (6)  we  quickly  arrive  at  the  result: 


F(eja)T)  = 


=  H(ejwT) 


(11) 


The  response  to  a  sampled  sinusoidal  input  can  be  readily  found  from  (11). 
e.  Representation  and  Geometric  Interpretation  of  System  Function 

We  see  that  the  system  function  H(z)  of  (8)  is  interpretable  as  a  frequency 
response  function  for  values  of  z  on  the  unit  circle  in  the  complex  z  plane.  Note 
that  the  radian  frequency  to  is  a  continuous  frequency  so  that  the  physical  signifi¬ 
cance  of  the  frequency  response  function  is  the  same  as  that  for  continuous  systems. 
Furthermore,  the  system  function  is  a  rational  fraction  whose  numerator  and  denom¬ 
inator  can  be  factored,  so  that  H(z)  is  uniquely  defined,  except  for  a  constant 
multiplier,  by  the  positions  of  its  poles  and  zeros  in  the  z  plane,  and  its  value 
for  any  point  z  determined  directly  by  the  distances  of  that  point  from  the  singu¬ 
larities.  We  are  thus  led  directly  to  the  geometric  interpretation  of  Fig.  5  whereby 
the  value  of  the  frequency  response  function  for  any  frequency  to  is  obtained  by 
rotating  by  an  angle  coT  about  the  circle,  measuring  the  distances  to  the  zeros 
Rj,  R2  •  •  •  • ,  the  distances  to  the  poles  P^,  P2  ....  and  then  forming  the  ratio,  so 
that  in  the  example  of  Fig.  5, 


(12) 


H(e^T)  algo  has  associated  with  it  a  phase,  which  is  given  by: 


0  -  0j  +  02  -  (0j  +  02  +  $3) 


(13) 


The  geometric  basis  for  digital  filter  design  is  thus  identical  in  principle 
with  the  geometric  basis  for  continuous  filter  design  with  the  following  single  and 
important  difference:  in  the  continuous  filter,  frequency  is  measured  along  the 
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imaginary  axis  in  the  complex  s  plane  whereas  in  the  digital  filter,  frequency 
is  measured  along  the  circumference  of  the  unit  circle  in  the  z  plane. 

f.  Several  Examples 

1st  order  difference  equation 
The  difference  equation: 

y(nT)  =  Ky(nT  -  T)  +  x(nT)  ( 14) 

has  the  system  function 

H(z)  =  - r  (15) 

1  -Kz 

The  frequency  response  function  is  then  obtained  by  letting  z  =  e^^  in 
(15),  which  yields, 


|  H(ejcoT)|  = 


J  1  +  K2  -  2K  cos  coT 


and  <p  =  angle  of  H(ejtoT)  =  -tan"* 1  cosS1^t  -  K 


(16) 


2nd  order  difference  equation 

The  difference  equation: 

y(nT)  =  Kx  y(nT  -  T)  +  K2  y(nT  -  2T)  +  x(nT) 

has  the  system  function 


(17) 


H(z)  = 


1 

IT  IT 

1  -Kx  z  -K2  z  z 


(18) 


and  the  resulting  frequency  response  magnitude 


|H(aJ“T)|  = 


T(1  —  cos  ooT  —  K2  cos  2coT)^  +  (K^  sin  toT  +  K2  sin  2ooT)^  J 


T/2 


(19) 
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and  phase 


Kj  sin  ooT  +  sin  2coT 


(20) 


1  —  K  ^  cos  coT  —  cos  2toT 


The  pole  positions  y  and  0  of  H(z)  are  given  by: 


(21) 


For  K2  negative  and  |  K2 1  >  —  the  poles  become  complex  conjugate  and 
as  they  approach  the  unit  circle,  the  well  known  resonance  effect  is  evidenced. 


-4 


Figure  6  is  a  plot  of  magnitude  vs.  coT  for  T  =  10 


a  resonant  frequency  of 


1500  cps  and  several  values  of  damping.  In  the  vicinity  of  the  poles,  the  magnitude 
function  behaves  very  much  like  a  continuous  resonator.  Of  course,  the  magnitude 
and  phase  must  be  periodic  in  coT;  this  follows  from  the  original  premise  of  a 
sampled  input  and  is  embodied  in  the  unit  circle  z  plane  representation.  We  state, 
without  proof,  the  well  known  facts  that  all  poles  of  H(z)  must  lie  within  the  unit 
circle  for  stability  and  that  poles  and  zeros  must  occur  on  the  real  axis  or  in 

th. 

conjugate  complex  pairs.  Also,  it  should  be  clear  from  equation  (8)  that  an  m 
order  difference  equation  has  a  z  transform  system  function  containing  m  poles 
and  r  zeros,  r  being  the  number  of  unit  delays  of  the  input. 
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SECTION  2 

a.  General  Discussion  of  Digital  Filter  Design  Techniques 

Since  much  information  is  available  on  continuous  filter  design,  a  useful 
approach  to  digital  filter  design  involves  finding  a  set  of  difference  equations 
having  a  system  function  H(z)  which  significantly  resembles  known  analog  system 

g 

functions.  The  work  of  Hurewicz  provides  a  technique  for  doing  this  in  an 
’impulse  invariant’  way.  By  this  we  mean  that  the  discrete  responses  to  an 
impulse  function  (see  section  Id)  of  the  derived  digital  filter  will  be  the  samples 
of  the  continuous  impulse  response  of  the  given  continuous  filter. 

Digital  filters  can  be  specified  from  a  desired  squared  magnitude  function 
using  procedures  akin  to  that  of  Butterworth  and  Chebyshev  continuous  filter  design 
procedure.  This  method  will  be  described  and  realizability  conditions  discussed 
in  an  Appendix. 

Another  technique  used  by  various  workers^’  ’  ’  uses  conformal 

mapping  to  transform  a  known  continuous  filter  function  into  a  digital  filter.  We 
will  refer  to  this  technique  as  ’frequency  invariant’,  since  critical  points  on  the 
continuous  filter's  magnitude  vs  frequency  curve  can  be  invariantly  transformed 
into  the  z  plane. 

Finally,  a  technique  referred  to  as  ’frequency  sampling’,  makes  use  of  the 

special  properties  of  an  elemental  digital  filter.  This  elemental  filter's  frequency 

sin.  x 

response  curve  resembles  a  — — —  function  and  has  a  linear  phase  vs  frequency 
characteristic.  By  suitably  combining  these  elemental  filters,  a  simple  design 
technique  for  a  large  variety  of  digital  filters  can  be  developed. 

Where  the  same  filtering  requirements  can  be  adequately  met  by  several 
different  digital  filters,  the  choice  between  them  depends  on  the  speed  of  execution 
of  a  computer  program  which  performs  the  difference  equation.  An  important 
factor  in  this  speed  is  the  number  of  multiplications  (not  counting  multiplications 
by  such  factors  as  ±  1  or  2n,  which  computers  can  do  quickly).  Some  digital  filters 
are  able  to  meet  essentially  the  same  requirements  as  others  with  substantially 
fewer  multiplications  per  output  sample,  and  these  are  to  be  preferred.  It  is 
important  to  stress  that  speed  of  execution  is  the  main  limiting  factor  in  the  utili¬ 
zation  of  digital  filtering  methods. 
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b.  Technique  1  -  Impulse  Invariance 


We  first  show  that  a  digital  filter  with  an  impulse  response  equal  to  the 
sampled  impulse  response  of  a  given  continuous  filter  can  be  derived  via  the 
correspondence 


m  A.  m  A. 

Y(s>  -  &  rrr  *  Si  — ^r-  =  H<2>  «2> 

1  i  1  J- 

1  —  e  z 

The  impulse  response  k(t)  of  a  continuous  filter  is  defined  as  the  inverse 
Laplace  Transform  of  its  system  function  Y(s),  given  in  general  form*  by  the  left 
side  of  (22).  Similarly,  the  impulse  response  h(nT)  of  a  digital  filter  is  defined 
as  the  inverse  z  transform  of  its  system  function  H(z),  which  can  be  expressed 
generally  by  the  right  side  of  (22). 

Thus 


i^m  A.  ^  m  -s.t 

k<'>  -  L  (s  1  sTs?)  -  Ai  e  <23> 

If  we  desire  that  h(nT)  =  k(t),  t  =  0,  T,  2T .  (24) 

m  -smT 

h(nT)  =  S1  A.  e  (25) 


*To  be  completely  general,  equation  (22)  should  also  include  terms  of  the  form 


A. 

l 

(s  +  si)'t/ 

which  correspond  to  terms  of  the  form 


f-1 


[fe? 


3<-l 


c 


A, 


1  -  e 


-aT 


s. 

i 
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Taking  the  z  transform  of  (25) 


H(z)  =  T  h(nT)  z 
n=o 


-n 


m  oo  -s.nT 

X  A .  X  e  z 
i=l  i  n=o 


m 


,  siT  -i 

1  —  e  z 


(26) 


Thus  the  condition  (24)  that  the  impulse  response  of  the  digital  filter  be 

equal  to  the  sampled  impulse  response  of  a  given  continuous  filter  Y(s)  leads  to  a 

digital  filter  defined  by  (26),  where  all  constants  A ^  and  s.  have  already  been 

specified  from  Y(s).  Thus,  by  means  of  the  correspondence  (22)  z  transforms 
13 

can  be  tabulated  . 

Example: 

The  simple  one  pole  RC  low  pass  filter  is  transformed  to  a  digital  filter  via 
the  correspondence: 


3.  v  cl 

s  +  a  ^  -aT  -T 

1  —  e  z 


(27) 


System  functions  of  various  resonant  circuits  may  be  expanded  by  partial  fractions, 
leading  to  the  correspondences 


s  +  a 

- 2 - 7 

(s  +  a)z  +  t> 


1  -aT  ,  ~  - 1 
1  —  e  cos  b  T  z 


,  0  -aT  ,  ~  -1  -2aT  -2 

1— 2e  cos  b  T  z  +e  z 


(28) 


(s  +  a)^  +  b^  1  —  2e  cos  b  T  z  *  +  e  z 


e  sin  b  T  z  1 

^T"  -2aT  -2 


(29) 


c.  Design  of  a  Digital  Lerner  Filter  which  is  "Impulse  Invariant" 
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The  Lerner  filter  is  defined  by  the  continuous  system  function 


m 

Y(s)  =  X 


B^s  +  a) 


with 


B1  =  7 


i=1  (s  +  a)  +  b^ 

B  -  (-I*™*1 
Bm  2 


B.  =  (— l)i+1  for  i  =  2, . . . ,  m  —  1 


(30) 
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and  the  pole  positions  (—  a  +  lh )  shown  in  Fig.  7.  It  has  been  shown  that  Lerner 
filters  have  a  high  degree  of  phase  linearity  and  reasonably  selective  pass  bands. 

From  the  correspondence  (28),  the  z  transform  of  the  digital  Lerner  filter 
is 


m 

H(z)  =  £ 


-flT  - 1 

B.(  1  —  e  ai  cosb.  Tz  ) 

i _ i _ ' _ 

■ii  o  "3-T  ,  r-p  “1  .  -2aT  -2 

1=1  1— 2e  cos  b.  T  z  +e  z 

l 


(31) 


Figure  8  shows  the  digital  realization  of  a  4th  order  band  pass  Lerner  filter. 
Each  of  the  four  parallel  sub  outputs  y.  is  computed  by  the  difference  equation 

y.(nT)  =  e'aT  cos  b.  T  [  2y.(nT  -  T)  -  x(nT  -  T)]  -  e"2aT  y^nT  -  2T)  +  x(nT)  (32) 


i  =  1,2, 3,4 


and  the  output  y(nT)  is  given  by: 


y(nT)  =  j  y1(nT)  -  y2(nT)  +  y3(nT)  -  j  y4(nT) 


(33) 


d.  Gain  of  Digital  Resonators 

The  correspondences  (28)  and  (29)  define  two  digital  resonators  which  are 
’impulse  invariants'  of  given  continuous  resonators.  In  practice,  digital  resonators 
can  be  specified  without  reference  to  continuous  resonators.  Specification  consists 
of  placing  the  pair  of  complex  conjugate  poles  and,  in  most  cases,  a  single  zero, 
from  which  the  difference  equation  can  be  quickly  derived.  Since  many  digital  filters 
are  simple  serial  or  parallel  combinations  of  these  resonators,  it  is  important  to 
understand  their  behavior. 

The  z  transform  of  a  resonator  with  poles  at  z  =  r  e^^r  and  a  zero  at  q 


is 


H(z)  =  - ■ - J-^T 

1  —  2r  cos  to  Tz  +  r  z 


(34) 


The  magnitude  vs  frequency  function  for  (34)  is  |H(e-*w^)|  and  can  be  written 
down,  by  inspection,  from  Fig.  5,  using  the  law  of  cosines. 
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|H(eJC0T)|  = 


1  +  q  —  2q  cos  COT 


[1  +  r2  -  2r  cos  (co  -  up)T]  [1  +  r2  —  2r  cos  (co  +  up)T] 


1/2 


(35) 


Case  1  q  =  cos  co  T 
-  r 

For  values  of  r  close  to  unity,  the  value  of  |H(e^^)|  at  the  resonance  co 
can  be  approximated  by 

lH<ei“T>l  -  iT-W  <36> 

which  is  independent  of  up.  Thus,  this  choice  of  q  makes  possible  the  design  of 
an  equal  gain  bank  of  resonators  (or  filters  composed  of  these  resonators)  covering 
a  wide  frequency  range. 

For  narrow  band  resonators  in  which  r  is  close  to  unity,  (36)  shows  clearly 
that  the  gain  at  resonance  is  usually  appreciably  greater  than  unity.  Knowledge  of 
filter  gains  is  required  for  determination  of  appropriate  register  word  lengths. 

Case  2  q  =  r  cos  co  T 

For  this  case,  the  difference  equation  for  (34)  can  be  written 

y(nT)  =  r  cos  C0rT(2y(nT  -  T)  -  x(nT  -  T))  +  r2  y(nT  -  2T)  +  x(nT)  (37) 

Execution  of  (37)  requires  only  2  multiplication  instructions  compared  to 
3  multiplications  for  general  q  and  this  case  is  thus  of  special  interest  for  real 
time  applications  and  when  total  computer  running  time  is  inordinately  lengthy. 
Sensitivity  of  the  resonant  gain  with  the  resonant  frequency  is  greater  than  for 
Case  1. 

Case  3  q  =  1 

This  case  yields  zero  gain  at  "d.  c.  ",  when  co  =  0,  which  is  often  desirable. 

j]~ 

As  in  Case  2,  two  multiplications  are  needed.  The  gain  for  upT  =  ps  JT  times 
the  gain  for  upT  close  to  zero. 

Case  4  q  =  0 

In  this  case,  the  zero  disappears.  In  the  design  of  formant  vocoders,  it  is 
desirable  to  design  resonators  without  zeros,  but  with  constant  gain  for  co  =  0, 
independent  of  up.  This  is  obtained  from  the  digital  system  function: 
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H(z)  = 


(38) 


2 

1  -  2r  cos  to  T  +  r 
r 

TI  2 

1  —  2r  cos  co  Tz  +  r  z 
r 

e.  Design  of  Digital  Filters  from  Continuous  Filters  Which  Have  Zeros  at 
Infinity 

A  large  class  of  analog  filters  are  defined  by  system  functions  of  the  form 


Y(s)  =  1  (39) 

n  (s  +  s  ) 

1=1  1 

where  the  denominator  is  a  product.  Such  filters  have  m  poles  at  finite  values  of 
th 

s  and  an  mu  order  zero  for  infinite  s.  Included  in  this  category  are  Butterworth, 
Chebyshev  and  Bessel  filters. 

In  order  to  design  impulse  invariant  digital  filters  based  on  (39),  the  pro¬ 
cedure  outlined  in  Section  2b  can  be  used;  Y(s)  is  expanded  in  partial  fractions, 
the  are  found  and  H(z)  is  obtained  from  the  correspondence  (22).  In  general, 
this  causes  zeros  to  appear  in  H(z),  although  there  were  no  finite  zeros  in  Y(s). 
However, ^hen  the  poles  at  — s^  are  close  to  the  imaginary  axis  in  the  s  plane,  such 
that  e  1  is  close  to  unity,  the  zeros  of  H(z)  can  be  ignored  and  H(z)  may  be 
approximated  by 


H(z) 


1 


m  x- 

xSi  0- 


-s.T 

l 

e  z 


(40) 


-sj 

In  fact,  it  can  be  shown  that  (39)  and  (40)  correspond  exactly  if  e 
can  be  replaced  by  1  —  s.T.  In  practice,  digital  band  pass  filters  several  hundred 
cycles  wide  of  the  Bessel,  Butterworth  or  Chebyshev  type  have  been  successfully 
programmed  for  10000  cycle  sampling  rates,  using  the  form  (40). 


Example  3  pole  Butterworth  Low  Pass  Filter 


The  system  function  of  the  continuous  filter  is 


Y(s) 


S1  s2  s3 

(s  +  S1)(s  +  S2XS  +  Sg) 


(41) 
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with 


S1  =  °c’  s2  =  K1+jyr>c’  S3  = 

tOc  being  the  cut-off  frequency,  defined  by  | Y ( j 60c )  |  =  .707.  Expansion  of  (41) 
into  partial  fractions  and  the  correspondence  (22)  leads  to  the  z  transform 


-«  T/2 


to 


H(z)  = 


-to  +  e 
c 


'  cos  C0CT  +  y|^  sin  C0cT  60cz  1 


-to  T 
1  c 

1  —  e  z 


_C0cT/2  1  ^  -1  'cVr 

1  —  2e  cos  4-/3  co  Tz  +e 


(42) 


which  has  the  diagrammatic  representation  shown  in  Fig.  9. 

If  a  cascade  representation  is  desired,  (42)  can  be  rewritten  to  give 


H(z)  = 


-1  -2 
C  z  +  Dz  z 


G- 


-to  T  N/' 

c  -_1\l-2e 


e  z  y 


-cocT 

“3T"  /3  ^  -1  -“cT  -2 

COS  — CO^T  Z  +  6  z 


0 


with  C 


“  “c[' 


-O)  T 
e  +  e 


-co  T 
c 


(yr  Sin4-“CT  -cos^l-  WCT 


•)] 


=  <oc[. 


-co  T 


-3co  T 
c 


and  D  =  co  I  e  —  e 


c-  -  2  - /I 


V3"  sin  -j-  «CT  +  cos  -  j- 


“cG] 


(43) 


H(z)  is  seen  to  have  two  zeros,  one  at  z  =  0  and  the  other  at  z  =  —  This 

zero  is  on  the  real  axis  and  increases  as  co  T  is  increased  from  zero.  Note 

c 

that  the  denominator  of  (43)  can  be  written  down  by  inspection  of  the  poles  of  the 

continuous  filter,  since  an  s  plane  pole  transforms  directly  into  a  z  plane  pole 
sT  D 

via  z  =  e  .  Thus,  the  zero  of  (43)  at  —  causes  the  only  error  in  the  assump¬ 
tion  that  an  s  plane  serial  system  transforms  directly  into  a  z  plane  serial  system. 

For  small  the  error  is  negligible,  since  the  zero  is  nearly  at  the  center  of 
the  z  plane  unit  circle  and  has  no  effect  on  the  magnitude  of  H(z).  Figure  10  shows 
a  plot  of  5-  vs  co  T.  We  see,  for  example,  that  the  zero  moves  about  5%  to  the 
right  of  the  origin  for  cocT  =  .1,  which  reduces  the  selectivity  of  the  cascade 
approximation  very  slightly. 
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f. 


Review  of  Butterworth  and  Chebyshev  Filters 


In  this  section,  we  review  briefly  design  procedures  for  continuous  Butterworth 
and  Chebyshev  filters,  which  will  be  needed  in  the  later  discussion  of  digital  filter 
design.  More  complete  treatments  are  available  in  standard  texts^’  ^ . 

The  Butterworth  filter  can  be  specified  by  the  relationship. 

|F(J0J)|2  =  - l- - j-  (44) 

1+(£) 

C 

where  coc  is  a  cut-off  frequency  and  |F(jui)|2  is  the  squared  magnitude  of  a  filter 

transfer  function.  The  poles  of  equation  (44)  lie  equally  spaced  on  a  circle  of  radius 

tOc  in  the  s  plane,  as  shown  in  Fig.  11.  For  n  odd,  there  will  be  poles  at  angles 

of  0  and  for  n  even,  the  first  pole  (beginning  at  angle  0)  occurs  at  an  angle 
7T 

of  2^.  It  can  be  shown  that  the  desired  transfer  function  F(s)  is  a  rational  func¬ 
tion  with  unity  numerator  and  denominator  determined  by  the  left  half  poles  of  Fig. 

11.  Plots  of  |F(jtd)|  of  (44)  for  several  values  of  n  are  shown  in  Fig.  12.  From 
these  plots  the  selectivity  properties  of  the  Butterworth  filter  become  clear. 

The  Chebyshev  filter  is  specified  by 


|F(jto)|2 


l  +  € 


(45) 


where  Vfl(x)  is  a  Chebyshev  polynomial  of  order  n  which  can  be  generated  by  the 
recursion  formula 

Vn+1(x)“2xVn(x)  +  Vn_1(x)  =  0  (46) 


with 


Vx(x)  = 


and 


V2(x)  = 


2x2  -  1 


The  Chebyshev  polynomial  has  the  property  of  equal  ripple  over  a  given 
range,  which,  with  added  specification  of  g,  leads  to  a  magnitude  function  of  a 
form  given  by  Fig.  13,  an  equal  ripple  in  the  pass  band  and  a  monotonic  decay  in 
the  stop  band.  The  ripple  amplitude  6  is  given  by 
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1 


(47) 


6  = 


The  poles  of  equation  (45)  lie  on  an  ellipse  which  can  be  determined  totally  by 
specifying  £,  n  and  cd,.  Figure  14  shows  this  ellipse,  the  vertical  and  horizontal 
apices  being  given  by:  btOc>  aco,  where 


b,  a 


■  m 


<E‘2  +  1  +  € 


l 

T\n 


7  ±  (J 


€"2  +  1  +  €' 


o'1} 


(48) 


where  the  b  is  given  for  the  plus  sign  and  a  for  the  minus  sign.  The  poles  on 
the  ellipse  may  be  geometrically  related  to  the  poles  of  two  Butterworth  circles  of 
radius  acoc  and  bcoc-  The  vertical  position  of  the  ellipse  pole  is  equal  to  the  ver¬ 
tical  position  of  the  pole  on  the  small  circle,  while  the  horizontal  position  is  that  of 
the  horizontal  position  of  the  large  circle  pole. 

g.  Technique  2  -  Digital  Filter  Specification  From  Squared  Magnitude  Function 

We  have  seen  in  Section  2f  that  the  Butterworth  and  Chebyshev  filters  are 
specified  by  choosing  suitably  selective  squared  magnitude  functions  such  as  in 
equations  (44)  and  (45).  The  same  procedure  is  possible  for  digital  filters  and  is 
described  in  this  section. 

Having  established  in  Section  Id  that  the  digital  filter  system  function  H(z) 
is  a  rational  fraction  in  z  \  it  follows  that  H(z)  for  z  on  the  unit  circle  is  a 
rational  fraction  of  e^^.  Thus,  the  squared  magnitude  |H(e^W^)|2  can  always  be 
expressed  as  the  ratio  of  two  trigonometric  functions  of  coT. 

An  example  of  a  squared  magnitude  function  suitable  for  low  pass  filtering  is 


|H(ejwT)|2  = 


tan 


1 

2n” 


1  + 


coT 

2 


tan 


2n 


co  T 

c 


(49) 


m 

Equation  (49)  is  plotted  in  Fig.  15  for  cOcT  =  ^  for  several  values  of  n. 
The  curves  obtained  are  similar  to  those  of  the  Butterworth  filter  of  equation  (44) 
plotted  in  Fig.  11.  The  cut-off  frequency  C0£  plays  the  same  role  in  both  con¬ 
tinuous  and  digital  case. 

Letting  z  =  e-*00^,  (49)  may  be  rewritten 
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(50) 


|H(z)|2  = 


tan 


o  w  T 
2n  c 


tan 


2n  “cT 


,  2n 

+(-‘)  (Itt) 


z  = 


We  see  that  (50)  is  a  rational  fraction  in  z  which  has  a  zero  of  order  2n  at 
—1.  The  poles  are  found  by  substituting,  in  (50) 


n  _  z  —  1 
P  z  +  1 


(51) 


from  which  we  can  ascg^un  that  the  2n  poles  of  p  are  uniformly  spaced  around 
a  circle  of  radius  tan  — ^ — .  The  ] 
formation  inverse  to  (51),  namely 


a  circle  of  radius  tan  — ^ — .  The  poles  of  z  are  then  readily  found  by  the  trans- 


z  =  1±P 
1  -P 


(52) 


Letting  p  =  x  +  jy  and  z  =  u  +  jv,  we  find  from  (52)  the  component 
relations 

,  2_  2  9 
u(x,  y)  =  - X  o  y~2  '•  v(x>  y)  =  - 9 - 7 

(i-xr  +  y  (i-xr  +  v 

The  circle  containing  the  poles  in  the  p  plane  satisfies  the  equation 


(53) 


2  2  2  ^0*^ 
xz  +  yz  =  tanz 


(54) 


From  (53)  and  (54)  it  can  be  shown  that  the  circle  of  (54)  maps  into  a  circle  in  the 
z  plane  centered  at  (u,v)  with  radius  p 


1  +  tan 


co  T 
2  c 


u  = 


— —  =  sec  co  T  ,  v  =  0 
co  T  c 


1  —  tan 


2  c 


2  tan 


co  T 

c 


P  = 


1  —  tan 


co  T 
2  c 


=  tan  co  T 
c 


(55) 
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For  odd  values  of  n,  the  2n  poles  of  p  have  the  x  and  y  coordinates 


co  T 

.  c  mTT 

x  =  tan  —n —  cos  - 

m  Z  n 


co  T 

_  c  mTT 

=  tan  -n —  sin  — — 
J  m  z  n 


m  =  0,  1 . 2n  —  1 


(56) 


For  even  values  of  n,  the  coordinates  are 


x  = 
m 


co  T 

.  c 1  2m  +  1 

tan  -2—  cos  — jjj —  it 


y_  =  tan 
J  m 


CO  T 
c 


m  =  0, 


1, .... 2n  —  1 


(57) 


cin 


2m 


From  (53)  and  (56)  the  corresponding  poles  in  the  z  plane  are  computed  to  be 


1  —  tan 


„  to  T 
2  c 


u  = 
m 


v  = 
m 


to  T  “  „  co  T 

i  „  .  c  mTT  2  c 

1—2  tan  —n —  cos  - +  tan  — n — 

z  n  z 


co  T 

o  .  c  mTT 

2  tan  —t, —  sin  - 

Z  n 


>  m  =  0,  1.  ..2n  — 1  (58) 


CO  T  ~  0  CO  T 

,  „  _  c  mtr  ,  „  z  c 

1—2  tan  — ft —  cos  - +  tan  — 

Z  n  Z 


Replacing  by  1  TT  yields  equivalent  formulas  for  n  even. 


Example  Find  the  poles  and  zeros  of  the  squared  magnitude  function  of  a 
low  pass  filter  with  3  db  attenuation  at  1250  cps  and  with  at  least  20  db  atten¬ 
uation  at  2000  cps.  Let  the  sampling  rate  be  10000  cps. 

The  cut-off  frequency  of  1250  cps  corresponds  to  tOcT  =  45°.  The  fre¬ 
quency  2000  cps  corresponds  to  coT  =  72  . 

The  squared  magnitude  function  (49)  becomes 
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1 


(59) 


|H(eJ‘oT))2  = 


1  + 


tan 


2n  coT 
2 


tan 


2n 


7T 

8 


The  appropriate  value  of  n  in  (59)  is  n  =  4,  obtained  by  setting 
COT  =  72°  and  |H(e^^)p  to  .01,  thus  satisfying  the  20  db  attenuation  condition. 
The  8  poles  in  the  p  plane  are  found  from  equation  (57).  Equation  (58)  can  now 
be  used  to  find  the  z  plane  poles  shown  in  Fig.  16,  which  also  shows  the  2n  zeros 
located  at  z  =  —1,  which  are  directly  derivable  from  equation  (50).  The  squared 
magnitude  function  is  thus  completely  specified  as  pole -zero  placements  in  the  z 
plane. 

In  the  Appendix,  an  analysis  is  made  of  the  necessary  relations  between  an 
assumed  squared  magnitude  function  and  the  digital  filter  specified  by  that  function. 

It  is  shown  that,  in  order  for  a  filter  to  be  realizable,  any  pole  inside  the  unit 
circle  (for  example  z^  in  Fig.  16),  must  have  a  mate  outside  the  unit  circle  of 
inverse  magnitude  and  the  same  angle.  Thus,  if  z^  =  r  e^,  there  must  be  a 
pole  (in  this  case  Zg)  given  by  ^-e^.  In  addition,  all  poles  must  occur  in  complex 
conjugate  pairs.  Therefore,  the  digital  filter  derived  from  Fig.  16  has  the  conju¬ 
gate  poles  z^,  Zg  and  Zg,  z^. 

The  above  argument  holds  for  the  zeros  as  well  including  the  special  case  of 
zeros  on  the  real  axis.  In  Fig.  16,  all  8  zeros  occur  at  z  =  —1.  The  derived 
filter  has  4  zeros  at  z  =  — 1. 

If  the  squared  magnitude  function  is  given  by 


|H(e>“T)|2  = 


(60) 


2  2 
1  +  *  Vn 


tan 

~r~ 

co  T 

tan 

c 

2  J 

1  — 

z  —  1 

f  _ 

z+  1 

lie  on  an  ellipse  in  the  p  plane 
which  has  the  same  properties  as  the  Chebyshev  ellipse  of  Fig.  14.  Using  the  nota¬ 
tion  of  Section  2f  and  Fig.  14,  the  p  plane  components  can  be  written 


WcT 

x  =  a  tan  — ^ —  cos  0 


co  T 

y  =  b  tan  — ^ —  sin  0 


(61a) 
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Substituting  (61a)  into  (53)  yields 


u  = 


!0- 


°°cT 

a  tan  — ^ —  cos  0  j 


o- 


CO  T 
c 

a  tan  — cos 


e)  +  b 


9  9  9 

tan  —75 —  sin  0 


-  1 


(61b) 


cocT 

2b  tan  — = —  sin  0 


v  = 


\  ,  u2  ”  2  WcT  .  2  ” 

(  1  —  a  tan  — ^ —  cos  0  j  +  b  tan  — ^ —  sin  0 


WcT 

Figure  17  shows  the  z  plane  mapping  for  a  tan  — ~ —  =  .  5  and 

b  tan  — ^ —  =  1.  The  ellipse  of  Fig.  14  maps  into  the  cardioid-like  curve  of 
Fig.  17  and  the  inner  circle  of  Fig.  14  maps  into  the  right-hand  circle  of  Fig.  17. 

The  outer  circle  of  Fig.  14  maps  into  a  circle  of  infinite  radius,  shown  by  the 
straight  line  of  Fig.  17.  The  points  shown  on  the  mapped  ellipse  of  Fig.  17  are 
the  actual  computed  points  from  equation  (61b). 

h.  Technique  3  -  Design  of  Digital  Filters  Using  Bilinear  Transformation  of 

Continuous  Filter  Function 

Another  approach,  using  design  directly  in  the  s  plane  will  now  be  considered. 
Suppose  we  have  a  stable  analog  filter  described  by  H(s).  Its  frequency  response 
is  found  by  evaluating  H(s)  at  points  on  the  imaginary  axis  of  the  s  plane.  If  in 
the  function  H(s),  s  is  replaced  by  a  rational  function  of  z  which  maps  the  imagin¬ 
ary  axis  of  the  s  plane  onto  the  unit  circle  of  the  z  plane,  then  the  resulting  H'(z), 
evaluated  along  the  unit  circle,  will  take  on  the  same  set  of  values  as  H(s)  evaluated 
along  the  imaginary  axis. 

Of  course,  this  does  not  mean  the  functions  are  the  same,  for  the  frequency 
scales  are  distorted  relative  to  one  another.  This  is  illustrated  by  the  simplest 
rational  function  which  maps  the  jto  axis  onto  the  unit  circle, 


s 


z  —  1 

Z  +  1 


(62) 


Let  the  analog  frequency  variable  be  co^,  and  let  the  digital  frequency  variable 
be  tdj-jT.  Then  the  functions  H(co^),  H'(tOj-jT)  take  on  the  same  values  for 
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oo^  =  tan 


ov  T 
D 


(63) 


Note  that  the  transformation  (62)  leads  to  a  ratio  of  polynomials  in  z.  Since 
it  maps  the  left  half  of  the  s  plane  onto  the  inside  of  the  unit  circle  we  can  be  sure 
that  it  will  always  yield  for  H*(z)  a  realizable,  stable  digital  filter. 

Equations  (62)  and  (63)  yield  a  technique  for  designing  a  digital  filter  by  analog 
techniques.  The  procedure  is  as  follows: 

(1)  Note  the  critical  frequencies  and  ranges  (pass  or  stop  band,  maximum 

attenuation  point,  etc. )  of  the  desired  digital  filter,  and  call  them  to^  T.  Compute 

a  new  set  of  frequencies,  to.  ,by 

Ai 


to.  =  tan 
A. 
i 


“Dj*1* 


(64) 


(2)  Design  a  transfer  function  H(s)  with  the  properties  of  the  digital 
filter  at  the  new  frequencies  and  ranges.  There  is,  of  course,  no  need  to  synthesize 
H(s). 

2  _  1 

(3)  Replace  s  by  .  .  in  H(s),  and  perform  the  algebra  necessary  to 
express  the  resulting  H'(z)  as  a  ratio  of  polynomials  -  this  yields  the  desired 
digital  filter. 

This  technique  is  illustrated  by  the  following  example,  in  which  we  design  a 

digital  filter  for  a  10  kc  sampling  rate,  which  is  flat  to  3  db  in  the  pass  band  of 

0  to  1000  cps,  and  which  is  more  than  10  db  down  at  frequencies  beyond  2000  cps. 

The  filter  must  be  monotonic  in  pass  and  stop  bands. 

From  Fig.  12  we  see  that  a  Butterworth  filter  meets  the  above  requirement  in 

the  analog  domain.  The  critical  frequencies  are  ton  T  =  2n  x  .  1  and  ton  T  = 

ul  u2 

2n  x  .2. 

(1)  Compute  co.  ,  to. 

A1  a2 


t oA  =  tan  2Tr  *  ‘  1  =  0.3249 
toA  =  tan  -  =  0.7265 
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(2)  We  design  a  Butterworth  filter  with  3  db  point  at  to  =  toA  = 

C  A  ^ 

.3249.  to  /to  =  2.235.  To  find  the  order  n  we  solve  1  +  (2.  235)2n  =  10 
a2  c 

and  obtain  n  =  2.  A  second  order  Butterworth  filter  with  0Oc  =  .  3249  has  poles 
at  s  =  .  3249  x  (—.  707  +  j  .  707)  =  — .  23  +  j.23,  and  no  zeros. 

H(S)  =  S1  ^2  _  2  x  .  232  _ .  1058 

(s  +  siHs  +  s2)  +  23)2  +  (.  23)2  s2  +  .  46  s  +  .  1058 

(3)  We  replace  s  by  p  yielding  H*(z). 


H*(z)  = 


.  1058 


CrTi~)  +  -46CItt)+  -1058 


_  .  1058(z  +  2z  +  1) _ 

1.  5658z2  -  1.  7884z  +  .  6458 


This  is  the  required  digital  filter.  It  requires  3  multiplications  per  output  point, 
or  if  DC  gain  can  be  tolerated,  only  two  multiplications. 

Problem:  Design  a  digital  filter  passing  from  0  to  100  cps  with  1/2  db  ripple, 

which  falls  off  monotonically  to  at  least  —19  db  at  183  cps.  Use  1000  cps  sampling 
rate. 

(1)  The  critical  frequencies  100  and  183  cps  are  transformed  to  analog 
frequencies. 


to 

c 


„  2  IT  X  100 

tan  T  x  1000 


tan  18°  =  .  32492 


to 

s 


tan 


2ir  x  183  _ 
2  x  1000 


.6498 


2t0 


c 


(2)  We  now  design  an  analog  filter  of  the  Chebyshev  variety.  One -half 

2 

db  ripple  corresponds  to  £  =  .  1220184.  To  find  the  required  order  we  solve 

2  2  19 

1  +  €  V  (to  /to  )  =  10  .  The  lowest  order  n  satisfying  this  relationship  is 

n  s  c 

n  =  3. 

A  unity  bandwidth,  1/2  db  ripple  Chebyshev  filter  is 
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constant 


Hl(s)  =  -7 


s  +  1.252913s  +  1.  5348954s  +  . 7156938 


Replacing  s  by  s/tdc  yields 


H(s)  = 


_ .0255842155 _ 

s+  .4127346s2  + . 16656307s  + . 0255842155 


where  the  constant  has  been  adjusted  for  unity  gain  at  s  =  0. 

z  *—  1 

(3)  We  replace  s  by  —  - -  . ,  giving  for  the  digital  filter  desired,  after 

2  *r  1  o 

multiplying  numerator  and  denominator  by  (z  +  1)  : 


it,/  \  _  .  02  55842  155  (z3  +  3z2  +  3z  +  1) 

^  \  o  r\  " 

1.  604882 18z  -3.  1694 18zz  +  2.44628634z  -  .  728243953 


which  is  the  required  digital  filter  design. 

It  is  worth  noting  a  useful  geometric  interpretation  of  Step  3.  Replacing  s 
2  _  1 

by  — -  .■  is  a  mapping  of  points  in  the  s  plane  onto  points  in  the  z  plane.  The 

Z  T  X 

following  short  table  gives  the  correspondence  of  some  critical  points  in  the  s 
and  z  planes 


s  plane 

z  plane 

0  + j  0 

1  + j  0 

00 

-1  + j  0 

0  +  j  1 

0  +  j  1 

0-j  1 

0-j  1 

-  1  + j  0 

0  +  j  0 

point  on  real  axis 

point  on 

point  on  imaginary  axis 

point  on 

point  on  any  line 

point  on 

A  very  similar  mapping  has  been  found  useful  in  some  applications,  and  graph 
paper  which  performs  the  transformation  can  be  purchased  under  the  name  "Smith 
Chart".  To  perform  the  mapping  of  our  application  we  take  a  conventional  Smith 
Chart  and  rotate  it  180°,  giving  a  chart  like  Fig.  18. 

The  location  of  any  point  -a  +  j  b  in  the  s  plane  (left  half)  is  used  to  find 
the  corresponding  location  in  the  z  plane  as  follows: 
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(1)  Locate  a,  a  point  along  the  centerline.  All  points  along  the  circle 
passing  through  this  point  correspond  to  points  in  the  s  plane  with  real  part  =  —a. 

(2)  Locate  b,  a  point  along  the  perimeter  of  the  outer  circle.  For  b  >  0 
use  the  top  semicircle;  for  b  <  0  use  the  lower  semicircle.  The  circular  arcs 
passing  through  this  point  correspond  to  points  in  the  s  plane  with  imaginary  part 
b. 

(3)  The  intersection  of  the  circle  with  the  circular  arc  is  the  point  in  the 
z  plane  corresponding  to  — a  +  j  b  in  the  s  plane. 

The  Smith  Chart  is  useful  when  the  function  in  the  s  plane  is  known  in  terms 
of  poles  and  zeros  or  poles  and  residues,  especially  when  the  z  plane  representation 
is  desired  in  that  form. 

Example:  Design  of  a  digital  high-pass  filter  with  stop  band  (0  -  500  cps) 
attenuation  greater  than  36  db  and  pass  band  (above  660  cps)  ripple  of  1.25  db. 
Sampling  rate  is  2.  5  kc. 

(1)  The  critical  frequencies  are  transformed  to  analog,  giving  as  the  stop 
band  limit,  .  72654  rad/sec,  and  as  the  pass-band  limit,  1.  09  rad/sec. 

(2)  The  specifications  are  met  with  a  fourth  order  elliptic  filter.  The 

17 

design  of  elliptic  filters  is  covered  in  the  literature  and  only  the  result  is  given 
here.  The  poles  are  at  s  =  —1.  1915812  +  jl.  5528835,  and  — .  11078815  + 
jl.  09445605,  and  the  zeros  are  at  s  =  0  (double),  and  +  jO.  69117051. 

These  are  located  on  a  Smith  Chart  in  Fig.  18.  Using  a  ruler  and  protractor 

+ i69  3° 

the  zeros  in  the  z  plane  are  found  to  be  at  z  =  1  (double),  and  e—  J  ’  and  the 
poles  are  found  to  be  at  z  =  .  586  e^  ^  ^  and  .  895  e—  ^  .  The  function  of 

z  defining  the  filter  is  thus 

H(z)  =  - fr-‘)V-7P*tH - 

(z  +  .  777z  +  .  3434) (zz  +  .  01877z  +  .  801) 

A  prime  difficulty  with  the  Smith  Chart  method,  illustrated  by  this  example, 
is  that  the  poles  cannot  be  located  with  great  accuracy.  In  compensation,  a  good 
deal  of  insight  into  the  operation  of  the  digital  filter  is  gained  in  the  course  of  the 
design.  Where  greater  accuracy  is  needed,  the  z  plane  poles  can  be  computed 
from  the  s  plane  poles  using 


z 


1  +  s 
1  —  s 
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which  is  the  inverse  of  equation  (62). 

2  _  1 

The  transformation  s  -» —  —  ^  is  not  the  only  rational  function  of  z  which 
maps  the  imaginary  axis  of  the  s  plane  onto  the  unit  circle.  For  example,  another 
such  transformation  is 


s 


(65) 


For  this  transformation  the  imaginary  axis  of  the  s  plane  maps  onto  the  top 

arc  of  the  unit  circle,  and  also  onto  the  bottom  arc.  The  origin  of  the  s  plane  maps 

,  .  +  i  ib  T 

onto  the  two  points  e—  J  ro  . 

Thus,  equation  (65)  maps  a  frequency  function  H(cd^)  into  H*(gOqT)  with 


co, 


cos  ipQ T  -  cos  tOpT 
sin  cOj-jT 


(66) 


This  means  that  equation  (65)  transforms  a  low-pass  H(to^)  into  a  band 
pass  H'(co^T).  Equation  (65)  can  thus  be  used  to  design  band-pass  digital  filters 
as  follows: 

(1)  Decide  on  a  center  frequency  i/)QT*  for  the  digital  filter.  This  may  be 
forced  by  the  specifications  or  it  may  be  available  to  simplify  the  choice  of  other 
parameters.  Compute  the  critical  frequencies  of  the  desired  analog  filter  from  the 
critical  frequencies  in  the  specification  of  the  digital  filter  using  Eq.  (66).  Equation 
(66)  will  often  yield  negative  frequencies,  which  is  alright  since  H(oo^)  will  be  an 

even  function  of  co.  . 

A 

(2)  Design  an  analog  filter  H(s)  with  the  translated  specifications.  It  is 
likely  that  one  or  more  of  the  specifications  will  be  superfluous. 

(3)  Replace  s  by 


s 


2 

z 


—  2z  cos  ipo T  +  1 


*A11  that  is  ever  really  needed  is  cos  i/)QT. 
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in  H(s),  and  perform  the  required  algebra  to  manipulate  the  resulting  H*(z)  into 
a  ratio  of  polynomials;  this  is  the  required  filter. 

One  example  will  probably  suffice  to  illustrate  the  method. 

Problem:  Design  a  digital  band -pass  filter  for  a  1000  cps  sampling  rate, 
to  pass  100  to  400  cps  with  ripple  free  attenuation  between  0  and  3  db.  At  45 
and  450  cps  the  filter  must  be  at  least  20  db  down,  and  must  fall  off  monotonically 
beyond  both  frequencies. 

(1)  *T  can  be  chosen  so  the  two  analog  3  db  points  are  negatives  of 
each  other.  If  the  digital  3  db  points  are  and  L^,  we  can  see  from  (66) 
that  if)Q T  should  be  chosen  so  that 

cos  i  (L,  +  L9) 

cos  ip  T  =  - i - - - —  (67) 

cos  j  (L1  -  L2) 

Using  (67)  yields  cos  =  0  for  this  case.  This  means  that  the  trans¬ 
formation  from  digital  to  analog  critical  frequencies  becomes 

60.  =  -cot  gclT 

A  D 

The  3  db  points  now  translate  to 

^db  =  ±  L3764 

The  20  db  points  are  not  equal  in  magnitude.  They  are 

60  =  -  3.  442 

S1 

60  =  3. 078 

s2 

(2)  The  problem  is  now  to  design  a  monotonic  filter  with  0-3  db  in  the 

region  0  <  60  <  1.  3764.  The  filter  must  be  20  db  down  by  60  =  3.  078,  which  will 

automatically  satisfy  the  requirement  at  60  =  3.  442.  A  Butterworth  design  seems 

to  be  called  for,  with  60  =  1.  3764. 

c 

For  —  =  2.  23  we  demand  20  db  attenuation.  We  calculate  the  order  n: 

60 

c 

1  +  (2.  23)2n  =  100 
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It  is  clear  that  n  =  3  will  suffice.  A  unity  bandwidth  Butterworth  filter  of  order 
3  is 


Hj(s) 


K _ 

77777777 


H(s)  is  found  by  replacing  s  by  T.  3764 


H(s)  = 


_  2.6075581 _ 

s3  +  2. 7528s2  +  3. 788954s  +  2. 6075581 


Z  -f  1 

(3)  To  find  the  required  digital  filter  we  let  s  =  — * -  in  H(s).  Note 

2  z  _  ]. 

that  since  this  is  going  to  yield  a  function  of  z  ,  the  resulting  digital  filter  will  be 
quite  simple. 


2. 6075581z6  -  7. 8226743z4  +  7. 8226743z2  -  2. 6075581 

10. 1493 12 lz6  -  5. 8588283z4  +  4. 2809203z2  -  0. 5714041 


If  the  gain  of  the  filter  is  not  important,  this  can  be  programmed  with  3 
multiplications  per  output  point,  plus  a  moderate  number  of  additions. 

It  should  be  obvious  how  the  techniques  of  this  section  can  also  be  used  to 
design  high-pass  and  band-elimination  filters, 
i.  Technique  4  -  Frequency  Sampling  Technique 

The  difference  equation 

y(nT)  =  x(nT)  —  x(nT  —  mT)  (68) 

has  the  z  transfer  function  1— z  m,  which  has  m  zeros  equally  spaced 
around  the  unit  circle,  at  points 

z ^  =  e-*2ir  m  ;  k  =  0,  1,  . . . . ,  m  —  1  .  (69) 

If,  in  equation  (68),  the  subtraction  is  replaced  by  an  addition,  the  transfer  func¬ 
tion  becomes  1  +  z  m,  for  which  the  zeros  are  also  equally  spaced  around  the 
unit  circle,  at 
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z,  =  e 


m 


k  =  0,  1, 


m 


-  1 


(70) 


The  magnitude  versus  frequency  curves  of  these  filters  repeat  with  period  — 
radians,  and  the  filters  are  often  referred  to  as  comb  filters.  These  filters  can 
be  incorporated  into  an  especially  appealing  type  of  design.  Before  presenting 
this  design,  let  us  examine  the  practical  significance  of  (68).  If  successive  values 
of  x(nT)  are  thought  of  as  contents  of  registers  in  computer  memory,  it  is  clear 
that  m  registers  must  be  set  aside  for  buffer  storage  in  order  to  execute  (68). 
This  is  usually  a  substantial  amount  of  memory  compared  to  that  needed  for 
execution  of  second  order  difference  equations.  Thus  possible  practical  use  of 
the  filters  we  are  about  to  describe  is  limited  to  systems  where  the  necessary 
memory  is  available,  or  to  systems  where  many  filters  share  a  common  input. 

It  is  important  to  note  that  the  delays  implied  by  (68)  can  be  effected  by  digital 
delay  lines,  which  are  presently  relatively  cheap  forms  of  memory. 

A  simple  resonator  can  be  placed  in  cascade  with  the  comb  filter.  Let  the 
resonator  be  of  the  type  considered  in  Section  2d,  case  4,  except  that  the  poles 
lie  directly  on  the  unit  circle.  Further  let  the  angle,  cor>  of  a  resonator  pole 
be  such  that  the  pole  is  coincident  with  a  zero  of  the  comb  filter. 


2ffk 

m 


cor  =< 


2tf(k  +  1/2)/  ^or  a  comk  filter  of  the 


m 


first 


second 


type 


(71) 


1  th 

The  poles  of  the  resonator  cancel  the  k  zero  of  the  comb  filter,  and  its  con¬ 
jugate.  In  what  follows,  we  shall  refer  to  the  resonator  which  is  used  to  cancel 
the  k1"*1  zero  as  the  k^1  elemental  filter. 

The  cascade  of  an  elemental  filter  with  a  comb  filter  is  a  composite  filter 
with  the  following  properties: 


(1)  The  impulse  response  is  of  finite  duration,  mT. 

(2)  The  magnitude  versus  frequency  response  is 


|H(e 


■  sin  mooT/2 
'cos  ooT  —  cos  co^T 


(72) 
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which  is  zero  at  all  the  radian  frequencies  for  which  the  comb  filter  is  zero, 

except  at  to  .  The  magnitude  at  to  is  5  esc  to  T. 

x*  r  £*  r 

(3)  The  phase  versus  frequency  is  exactly  linear  except  for  discontinu¬ 
ities  of  it  radians.  These  discontinuities  occur  where  the  magnitude  response  is 
zero. 

(4)  The  phase  difference  between  two  composite  filters  with  resonant 

frequencies  oo  and  is  V  for  tor  <  to  <  and  is  zero  outside  these 

bounds. 

(5)  The  amplitude  of  any  composite  filter  is  zero  at  the  resonant  fre¬ 
quencies  of  all  the  other  composite  filters. 

As  m  is  made  large,  the  magnitude  response  of  a  cascaded  filter  becomes 

like 


P  sin(tO  -  tOr)T 
L  (tO-tOr)T 


sin(to  +  tor)T 

(to  +  topT  J 


in  shape.  These  properties  suggests  that  any  desired  magnitude  response  could  be 
obtained  by  adding  together  the  weighted  outputs  of  cascaded  comb  and  elemental 
filters,  just  as  any  "band-limited"  time  function  can  be  formed  from  a  weighted  sum 
of  delayed  functions.  Let  us  examine  this  idea,  which  we  call  frequency 

sampling,  in  some  detail. 

A  sufficiently  "narrow  band"  frequency  response  function  (one  for  which  the 
frequency  response  is  a  sufficiently  smooth  function  of  frequency)  is  sampled  at  equally 
spaced  points,  with  radian  frequencies: 


u  =  2nk  or  2g(k  +  1/2) 

km  m 


k  =  0,  1,...,  m  —  1 


depending  on  which  kind  of  comb  filter  will  be  used.  Let  the  sample  value  of  the 
amplitude  at  frequency  tOj,T  be  W^.  An  elemental  filter  of  resonant  frequency 
to^,  cascaded  with  a  comb  filter  of  delay  mT  and  a  gain  of  W^.  sin  co^T  are  used 
to  provide  an  "elemental  frequency  response"  of  Wj,  at  radian  frequency  oo^T  and 
zero  at  the  other  sampling  frequencies. 

Since  the  phases  at  resonance  of  the  consecutive  elemental  filters  differ  by  17, 
the  gains  of  the  odd  numbered  elemental  filters  are  to  be  multiplied  by  —  1.  The 
desired  input  to  the  filter  is  applied  to  the  comb  filter,  which  is  shared  among  all 
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the  elemental  filters,  followed  by  the  gains  (and  sign  changers  for  odd  numbered 
elemental  filters).  The  outputs  of  all  the  elemental  filters,  with  proper  gains,  are 
added  together  to  give  the  desired  filter  output.  The  resulting  filter  has  an  impulse 
response  of  duration  mT,  a  frequency  response  with  linear  phase,  and  an  amplitude 
response  which  agrees  with  specifications  at  the  sampling  frequencies  and  connects 
the  sampling  points  smoothly.  The  variety  of  filters  which  can  be  programmed  with 
this  technique  is  quite  large. 

There  are  some  practical  problems  to  be  considered  before  the  frequency 
sampling  method  is  applied.  For  one  thing,  the  resonant  poles  of  an  elemental  filter 
cannot  exactly  cancel  the  zeros  of  a  comb  filter  because  of  quantization.  Thus  it  is 
wise  to  move  both  the  zeros  of  the  comb  filter  and  the  poles  of  the  elemental  filters 

T  2^ 

slightly  inside  the  unit  circle,  with  a  radius  of  something  like  e  =1—2 

We  have  successfully  programmed  filters  with  the  poles  and  zeros  at  radii 
-12  -27 

ranging  from  1—2  to  1  —  2  ,  with  little  change  in  the  behavior  of  the  filter. 

Another  change,  while  not  necessary,  is  useful  to  keep  in  mind  for  band-pass 

filters.  In  the  pass  band  it  is  common  for  the  samples  to  be  equal.  Thus  it 

would  be  desirable  if  all  the  elemental  filters  had  the  same  gains  at  resonance.  If 

-aT 

a  zero  is  put  at  (cos  td^T)  e  ,  the  gain  of  each  elemental  in  cascade  with  the 
comb  filter  becomes  m/2.  This  has  a  slight  effect  on  the  magnitude  and  phase 
response,  negligible  for  large  m.  The  modified  comb  filter  z  transform  thus 
becomes 


H(z)  = 

and  the  modified  kcn 


,  -maT  -m 
1  —  e  z 

elemental  filter  becomes 


Hk<2) 


-aT  - 1 

1  -  e  cos  tOj/T  z 

~i  o  -  aT  - 1  -2aT  -  2 

1  —  2  e  cos  C0j/T  z  +  e  z 


(73) 


(74) 


It  is  worth  noting  that  the  introduction  of  the  additional  zero  does  not  require 
another  multiplication  since  twice  the  numerator  coefficient  is  also  present  in  the 
denominator.  The  response  of  a  filter  of  the  form  of  equations  (73),  (74)  is  shown 
in  Fig.  19. 


31 


Example:  Bank  of  Band  Pass  Filters 


It  was  desired  to  design  a  bank  of  band-pass  filters  (with  a  common  input), 
each  400  cps  wide,  covering  the  band  300  cps  to  3100  cps.  The  filters  are  to  be 
as  selective  as  possible  but  with  minimum  ringing  time.  A  further  requirement  is 
that  the  contiguous  filters  cross  at  — 3  db  of  the  midband  gain.  None  of  the  standard 
designs  had  satisfactory  selectivity  combined  with  short  ringing. 

The  filters  chosen  were  frequency  sampling  filters,  each  composed  of  seven 
elemental  filters.  The  consecutive  zeros  were  100  cps  apart.  Since  the  sampling 
rate  was  12.5  kc,  m  was  125.  The  general  form  of  such  a  filter  has  z  transform 


H(z)  =  Ql-e 


r+6 

-maT  -m'NV 

2  ^ 
k=r 


(-if 


.  -aT  2irk  -1 

1  —  e  cos - z 

m 

0  -aT  2irk  -1  -2aT  ; 

2e  cos - z  +  e  z 

m 


2  Wk 


(75) 


The  design  of  the  filter  is  completely  specified  by  choosing  r  and  the  set  of 
in  (75).  Since  3  db  crossovers  400  cps  apart  were  required,  and 

were  chosen  to  be  .707.  The  three  center  terms  had  W^,  =  1.  The  end 
term  gains,  Wr  and  were  found  empirically  to  be  .221  for  satisfactory  out 

of  band  rejection. 

The  design  of  the  300-700  cps  filter  is  illustrated  in  Fig.  20,  and  the  exper¬ 
imental  frequency  response  of  it  and  the  next  higher  filter  (700-1100  cps)  are  shown 
in  Fig.  21. 
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Quantization  Noise  of  Digital  Filters 

When  a  digital  filter  is  realized  with  a  digital  arithmetic  element,  as  is  the 
case  with  a  computer,  additional  considerations  are  necessary  to  describe  the  per¬ 
formance  of  the  filter.  There  are  three  obvious  degradations:  1)  quantization  of 
the  coefficients  of  the  difference  equation,  2)  quantization  of  the  input,  and  3)  quan¬ 
tization  of  the  results  of  computations.  The  first  two  effects,  although  important, 
are  simple  to  understand.  Quantization  of  the  coefficients  changes  them  to  slightly 
different  coefficients.  This  happens  once  and  for  all,  resulting  in  a  new,  slightly 
different  filter.  For  very  high  Q  filters  this  may  be  important;  even  instability 

may  result.  Quantization  of  the  input  is  equivalent  to  adding  a  noiselike  term  which 

18 

has  been  well  described  in  the  literature.  We  rely  heavily  on  the  results  of  Bennett 

in  treating  this  and  similar  effects  of  quantization  in  a  statistical  manner. 

Quantization  of  results  of  computation  is  more  complicated  because  these 

results  are  used  in  later  computations,  and  thus  the  effect  of  the  error  on  future 

computations  must  be  understood.  Let  us  first  note  that  quantization  of  the  results 

cannot  be  avoided  because  after  each  iteration  of  the  difference  equation,  the  number 

of  bits  required  for  exact  representation  of  the  result  increases  by  about  as  many 

bits  as  in  the  representation  of  the  coefficients.  Thus,  without  quantization,  the 

number  of  bits  required  would  grow  without  bound. 

In  the  following,  we  define  quantization  as  the  replacement  of  the  exact  value 

of  a  quantity  by  the  value  of  the  nearest  of  a  set  of  levels  differing  by  steps  of  EQ. 

Thus,  the  maximum  error  introduced  by  a  quantization  is  E  /2,  and,  except  for  a 

2  0 

few  degenerate  cases,  the  mean  squared  error  is  Eq/12. 

Figure  22  shows  a  representation  of  a  digital  filter  in  which  the  output  is 
quantized  before  being  used  in  further  iterations. 

We  can  write  a  "linear”  difference  equation  describing  this  device  if  we 
represent  the  quantization  error  by  v(nT). 

N  M 

r(nT)  =  v(nT)  +  £  a.  x(nT  -  iT)  +  £  b.  r(nT  -  iT)  (76) 

i=o  i  i=l  i 

This  equation  has  a  z  transform  form: 
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(77) 


N  . ,  M  • 

R(z)  =  V(z)+X  a.z  X(z)  +  S,  b.  z  R(z) 
1=0  i  i=l  i 


(78) 


We  notice  that  the  z  transform  of  the  output  consists  of  two  terms,  one  of 
which  is  the  desired  output.  The  other  term  depends  only  on  the  poles  of  the  filter 
and  on  the  noise  sequence  v(nT). 

In  what  follows,  assume  that  the  input  to  the  filter  is  of  such  level  that  the 

output  with  no  quantization  noise  has  a  mean  squared  value  of  unity.  Assume  that 

this  output  is  large  compared  to  E  ,  so  that  the  samples  v(nT)  are  independent 

0  18 

and  uniformly  distributed  from  —  EQ/2  to  Eq/2  .  Let  us  call  the  denominator 
of  the  right  side  of  (78)  Q(z).  Our  aim  is  to  determine  the  mean  squared  value  of 
the  noise  term  of  (78)  in  terms  of  Eq,  the  quantization  step  size. 

13 

For  this  determination,  the  two-sided  z  transform  of  the  autocorrelation 
function  of  a  sequence  is  needed.  If  a  one-sided  sequence  has  z  transform  A(z) 
then  its  autocorrelation  function  has  a  two-sided  z  transform  A(z)  A(— ).  It 
follows  that  a  sequence  with  z  transform  V(z)/Q(z)  has  an  autocorrelation  func¬ 
tion  whose  z  transform  is  the  ratio  of  the  z  transforms  of  the  autocorrelation 
functions  of  v(nT)  and  q(nT). 

We  are  really  interested  in  a  normalized  autocorrelation  function, 


which,  for  v(nT)  is  the  sequence 


Ay(nT)  =  0,  ...,0,  E^/12,  0 
2 

with  z  transform  A  (z)  =  E  /12. 

vv  '  o' 


0 


For  the  noise  output,  the  normalized  autocorrelation  function  has  z  transform: 
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The  mean  squared  output  noise  is,  of  course,  given  by  the  coefficient  of  z° 
in  the  z  transform  of  its  autocorrelation  function.  For  a  two-sided  z  transform 
this  coefficient  is  given  by 


M.  S.  O.  N.  =  —£5-  (  )  residues  of  - - - (79) 

^  ^  zQ(z)  Q(~) 

where  only  the  residues  at  poles  inside  the  unit  circle  are  to  be  included  in  the  sum. 
Example: 

Consider  a  one  pole  filter  with  H(z)  =  - - — T.  Then  — w  L, .  ■  1  has 

1-  az  z  Q(z)  Q(l/z) 

a  pole  at  z  =  a  and  a  residue  of  — ^  at  that  pole.  The  mean  squared  output 

noise  for  a  digital  realization  of  this  filler  is  proportional  to  this  residue,  which  is 

plotted  versus  pole  position  in  Fig.  23,  with  the  constant  of  proportionality  being 

E  2/12. 
o 

Two  filters ,  one  using  36  bit  arithmetic  and  one  using  a  variable  number  of 
bits  fewer  than  36,  were  simulated,  with  a  common  input  to  both  filters,  and 
with  H(z)  as  above.  The  36  bit  case  was  considered  unquantized  and  the  differ¬ 
ence  in  the  filter  outputs  was  thus  a  good  approximation  to  the  quantization  noise 
introduced  by  the  less  precise  filter.  The  mean  squared  value  of  this  noise  was 
measured  for  a  random  noise  input  with  28  and  29  bit  arithmetic,  and  for  a 

sine  wave  input  at  about  0.  1  of  the  sampling  frequency,  using  29  bit  arithmetic. 

2 

The  measured  results,  normalized  with  respect  to  Eq/12,  are  shown  along  with 
the  theoretical  result  in  Fig.  23. 

We  have  arrived  at  a  very  fortunate  expression  since  we  can  quickly  solve  it 
for  (— log2  Eq)  which  is  the  number  of  bits  of  quantization  needed  below  the  output 
signal  level,  for  a  given  amount  of  mean  squared  output  noise.  Say  we  require 
that  the  mean  squared  output  noise  be  F  db  below  a  unit  output  signal  level.  Then 

E2 

F  =  10  log1Q  (r  residues )j-  (80) 


b  ,  =  —  log9  E  =—.  166  F  —  1.  79  +  1.  66  log. n(  S residues  of - 7-  )  (81) 

_i  1  °  UV  zQ(z)Q(i)y 

b  ^  is  the  number  of  bits  which  must  be  retained  below  the  unit  signal 
level;  an  additional  3-5  bits  ought  to  be  retained  above  the  unit  signal  level  to  protect 
against  overflow. 
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The  preceding  description  of  quantization  noise  is  valid  for  digital  filters 
realized  in  the  form  of  Fig.  1.  If  the  realization  is  in  the  canonical  form  of  Fig.  4 
then  quantization  noise  is  added  to  w(nT)  rather  than  to  r(nT),  since  it  is 
w(nT)  which  is  used  in  further  computations.  Figure  24  shows  the  digital  arrange¬ 
ment  of  the  canonical  form.  We  write  the  two  difference  equations: 

M 

w(nT)  =  x(nT)  +  £  b.  w(nT  —  iT)  +  v(nT)  (82) 

i=l  1 


r(nT)  =  §  a.  w(nT  -  iT) 
i=o  i 

These  may  be  z  transformed,  giving 

W(z)  =  (X(z)  +  V(z))/Q(z) 

R(z)  =  P(z)  W(z) 

N  _i 

where  P(z)  =  £  a.  z 

i=o  1 


(83) 


(84) 

(85) 


Combining  these  into  a  single  equation 

R(Z>  -  X(z)  +  V(z)  gg-  (86) 

we  see  that  in  the  canonical  representation  the  quantization  noise  is  filtered  by 
both  the  poles  and  the  zeros  of  the  filter.  If  desired,  the  procedures  leading  to 
equations  (79),  (80)  and  (81)  could  be  repeated,  giving  as  the  mean  squared  output 
noise 


E  2  P(z)  P(-) 

M. S. O. N.  =  £  residues  of  Q - ^~Tj 


'zQ(z)  Q £) 


(87) 


and  giving  the  number  of  bits  required  below  a  unit  signal  level  as 


P(z)P0  . 

b^  =—  .166F  —  1.  79  +  1.66  log-^QC  residues  of - - - p  ) 


zQ(z)  Q0 


(88) 
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APPENDIX 


REALIZABILITY  OF  DIGITAL  FILTERS  WITH  CERTAIN  SQUARED 

MAGNITUDE  FUNCTIONS 

In  the  following  we  shall  find  it  useful  to  define  polynomials  of  the  form 
k  .  k-1  , 

az  +  a,z  +  ...+a,z  +  a„  , 
o  1  1  o 

r  k  ~  r 

such  that  the  coefficient  of  z  is  equal  to  the  coefficient  of  z  ,  as  mirror 
image  polynomials  of  order  k,  and  to  summarize  some  of  the  properties  of  the 
polynomials. 

(1)  A  mirror  image  polynomial  of  order  k  has  roots  which  occur  in 
reciprocal  pairs  for  k  even. 

Proof:  Let  z  be  replaced  by  1/z  in  the  original  M.  I.  P.  and 
equate  it  to  zero. 

-k  -k+1  -1  „ 

az  +a,z  +a,z  +a  =  0 
o  1  1  o 

If  we  multiply  by  z  ,  which  introduces  no  new  roots  to  the  equation,  we  see  that 
the  roots  of  the  polynomial  in  1/z  are  the  same  as  the  roots  of  the  polynomial  in 
z.  Thus,  the  roots  either  occur  in  reciprocal  pairs,  or  are  self- reciprocals.  For 
k  odd,  one  of  the  roots  cannot  be  part  of  a  reciprocal  pair.  This  root  is  its  own 
reciprocal,  and  since  it  must  be  real,  it  is  either  +1  or  —1.  We  will  generally 
not  be  interested  in  odd  order  M.  I.  P.  *s. 

(2)  The  sum  of  mirror  image  polynomials  of  the  same  order  is  a  M.  I.  P. 
of  that  order. 

(3)  A  mirror  image  polynomial  of  order  k  plus  z  times  a  mirror 
image  polynomial  of  order  (k  —  2r)  is  a  mirror  image  polynomial  of  order  k. 

(4)  The  product  of  mirror  image  polynomials  is  a  mirror  image  poly¬ 
nomial. 

(5)  A  polynomial  whose  roots  occur  in  reciprocal  pairs  is  an  M.  I.  P.  (5) 
is  proved  by  multiplying  together  the  factors  which  are  reciprocals,  noting  that 
each  product  is  =  M.  I.  P. ,  and  applying  (4). 

We  shall  use  the  above  properties  of  mirror  image  polynomials  to  prove  the 
existence  of  digital  filters  with  certain  squared  magnitude  functions.  Suppose  a 
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squared  magnitude  function,  F(coT),  is  given.  We  can  always  replace  coT  by 
—  j  log  z,  equivalent  to  z  =  to  give  a  function  G(z)  which  is  equal  to 

F(toT)  when  evaluated  along  the  unit  circle.  For  our  purposes  G(z)  must  be  a 
rational  function  of  z  with  real  coefficients.  We  can  guarantee  this  if,  in  the 
squared  magnitude  function,  coT  only  appears  as  part  of  the  following  expressions: 

.  2  toT  2  toT  2  COT  .2  COT  2  2  toT  ,  . , 

sin  -y-,  cos  ~2~>  ^n  cot  ,  sec  -j-,  esc  — ,  and  if  these  expres¬ 

sions  are  combined  only  by  addition,  subtraction,  multiplication  and  division.  This 

coT 

is  because  the  squared  triginometric  functions  of  —  yield  rational  functions  of 

z,  and  rational  functions  form  a  field  with  multiplication  and  addition.  The  corres- 

coT 

pondence  of  squared  trigonometric  functions  of  —  with  rational  functions  of  z 
is  given  below: 


TABLE  1 


sin 


2  coT  (z  -  1)^ 
2  "  — 4z 


cos 


2  toT  (z  +  1)' 
2  4z 


tan 


2  coT  -(z  -  1) 

-  “* - T 

(z  +  l)2 


cot 


2  coT  _  — (z  +  1)" 

2  ,  IX2 
(z  -  1) 


sec 


2  coT 


4z 


2  (2  +  l)2 


CSC 


2  toT 


-4z 


(z  - 1  y 


Now  let  us  assume  that  we  have  been  given  an  F(toT)  with  the  preceding 
restrictions,  found  G(z),  and  found  all  the  poles  and  zeros  of  G(z).  To  find  a 
transfer  function,  H(z),  with  squared  magnitude  function  F(toT),  we  should  like 
to  use  the  following  procedure: 

(1)  Discard  any  poles  or  zeros  at  the  origin. 
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(2)  Put  the  remaining  poles  in  one  to  one  correspondence  with  each  other 
in  such  a  way  that  for  each  pair  of  poles  the  distance  from  one  pole  to  any  point  on 
the  unit  circle  is  in  constant  ratio  to  the  distance  from  that  point  on  the  unit  circle 
to  the  other  pole  of  the  pair. 

(3)  Repeat  step  2  for  the  zeros. 

(4)  Discard  one  of  each  corresponding  pair  of  poles  and  zeros ,  making 
sure  that  the  retained  poles  are  all  within  the  unit  circle. 

(5)  Form  H(z)  from  the  retained  poles  and  zeros  only. 

This  procedure  works  only  if  steps  2  and  3  are  possible.  The  magnitude 
of  H(z)  is  proportional  to  the  product  of  the  distances  from  the  unit  circle  to  retained 
zeros  divided  by  the  product  of  distances  to  the  retained  poles,  and  thus  is  propor¬ 
tional  to  the  square  root  of  F(toT).  Note  that  the  singularities  of  G(z)  at  the  origin 
contribute  nothing  to  a  magnitude  function  since  the  distance  from  the  origin  to  the 
unit  circle  in  always  one. 

Steps  2  and  3  are  possible  if  the  poles  (and  zeros)  occur  in  pairs  whose 
locations  are  equal  or,  more  commonly,  conjugate  reciprocals  of  each  other, 

i  £  |  if) 

reJ  and  —  eJ  .  Since  the  singularities  of  G(z)  occur  in  conjugate  pairs  it  is 

enough  to  require  the  roots  to  occur  in  equal  or  reciprocal  pairs.  By  the  following 

proof  we  show  that  the  ratio  of  distances  from  any  point  on  the  unit  circle  to  two 

singularities  which  are  conjugate  reciprocals  is  constant.  Proof:  Let  the  two  fixed 

singularities  be  reJ  and  —  e  .  By  the  theorem  of  Appollonius  the  locus  of  points, 

the  ratio  of  whose  distance  to  two  fixed  points  is  constant,  is  a  circle  whose  center 

id 

is  colinear  with  the  two  fixed  points.  Consider  the  ratios  at  the  two  points  eJ 
and  The  ratios  are  equal  to  r  at  both  points,  and  the  only  circle  passing 

through  both  points  whose  center  is  colinear  with  the  two  fixed  points  is  the  unit 
circle.  Thus,  the  unit  circle  is  the  circle  of  Apollonius,  and  the  ratio  of  the  dis¬ 
tance  to  the  singularities  is  constant. 

The  corresponding  proof  for  equal  singularities  is  trivial. 

If  after  discarding  poles  or  zeros  at  the  origin  the  numerator  and  denominator 
of  G(z)  become  mirror  image  polynomials,  it  is  possible  by  following  the  five  step 
procedure  given  above  to  find  H(z)  with  squared  magnitude  function  F(cOT)*.  We 

*This  is  not  strictly  correct  if  the  mirror  image  polynomials  have  simple  roots  on 
the  unit  circle,  which  are  then  their  own  conjugate  reciprocals.  Subsequent  state¬ 
ments  should  be  qualified  by  taking  this  possibility  into  account.  However,  this  is 
a  limiting  case,  unlikely  to  occur. 
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shall  now  enumerate  a  few  types  of  F(cdT)  which  yield  mirror  image  polynomials 
for  numerator  and  denominator  of  G(z). 

2 

Consider  F(OOT)  a  rational  function  of  sin  -j— .  Then  G(z),  by  substitu¬ 
tion  from  Table  1,  is  of  the  form: 


G(z) 


(2-1) 

— 4z 


(2~1) 

— 4z 


2 


2 


1 

1 


n 


m 


♦  4^] 


2-,n- 1 


+ 


2  _m- 1 


+  . . . . 


which  can  be  put  in  the  form 


G(z) 


(aQ(z  -  l)* 2)n  +  aA((z  -  l)2)n' \  —4z)  +  a2((z  -  l)2)n_2(-4z)2  +  . . .  )(-4z)m'n 
b0((z  -  l)2)m  +  bx((z  -  l)2)m_1(-4z)  +  b2((z  -  l)2)m  2(— 4z)2  +  . . . 


The  (—  4z)  contributes  only  poles  or  zeros  at  the  origin  and  is  discarded. 

2 

Note  that  powers  of  (z  —  1)  are  mirror  image  polynomials  by  property  (5),  and 
thus  the  entire  numerator  of  G(z)  is  a  mirror  image  polynomial  by  property  (3). 
The  same  reasoning  holds  for  the  denominator. 

2  coT 

Exactly  similar  reasoning  shows  that  rational  functions  of  cos  also  yield 

2  coT  ^ 

realizeable  digital  filters.  Rational  functions  of  tan  —  yield  G(z)  of  the  form 


G(z) 


aQ(z  —  l)2n  +  aj(z  —  l)2^n  ^(z  +  l)2  +  . . .  +  an(z  +  l)2n 

bQ(z  -  l)2m  +  bj(z  -  l)2(m'  ^(z  +  l)2  +  . . .  +  bm(z  +  l)2m 


(z  +  l)2(m-n) 


Here  we  note  by  repeated  application  of  property  (4)  that  again  both  numerator 
and  denominator  are  mirror  image  polynomials  and  therefore  a  digital  filter  exists 


with  the  desired  squared  magnitude  function.  Again  exactly  similar  reasoning 

t2  ooT  2  ooT  , 

-x — ,  sec  -x —  and 

2  ouT 

CSC  — R— . 


suffices  to  extend  the  proof  to  rational  functions  of  cot'* 


Similar  discussions  show  that  rational  functions  of  sums  of  products  or  ratios 
of  the  squares  of  the  trig  functions  of  Table  1,  i.e. ,  (sin2  ^r-+  tan^  ^-)  also 
yield  G(z)  which  is  a  ratio  of  mirror  image  polynomials.  Further  extensions  can 
be  made,  but  we  shall  be  content  here  with  these  few  examples  and  techniques. 
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Fig.  1.  Pictorial  representation  of  m  order  linear  difference  equation. 
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Fig.  2.  Serial  arrangement  of  difference  equations. 
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Fig.  3.  Parallel  arrangement  of  difference  equations. 


Alternate  representation  of  digital  filter  (drawn  for  m  =  r) 


Fig.  5.  Geometric  interpretation  of  digital  filter  frequency  response 
function. 
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Fig.  6.  Response  of  digital  resonator  for  several  values  of  damping. 
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Fig.  7.  Pole  positions  for  Lerner  filter. 
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Fig.  8.  Configuration  of  4LU  order  (8  pole)  digital  Lerner  filter  (dotted 
box  shows  configuration  of  pole  pair  G.(z)). 
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Fig.  9.  Impulse  invariant  3 -pole  Butterworth  low  pass  filter. 
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Fig.  10.  Graph  of  showing  position  of  extra  zero  introduced  by  impulse 
invariant  design  technique. 
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Fig.  11.  Foies  of  Butterworth  squared  magnitude  function,  equation  (44)  for 
n  =  5.  Poles  are  equally  spaced  with  angular  separation  36°.  The  system 
function  Y(s)  is  defined  by  the  poles  in  left  half  plane. 
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Fig.  12.  Magnitude  vs  frequency  of  Butterworth  filters. 
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Fig.  13. 


Magnitude  vs  frequency  of 


4th 


order  Chebyshev  low  pass  filter. 
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Fig.  15.  Magnitude  vs  frequency  of  frequency  invariant  digitized  Butterworth 
filters. 
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Fig.  16.  Poles  and  zeros  of  squared  magnitude  function  for  example  of 
section  2g. 


Fig.  17.  Bilinear  mapping  of  Chebyshev  ellipse  into  z  plane  for 

to  T  to  T 

c  o 

a  tan  — ^ —  =  .  5,  b  tan  — ^ — 

degree  increments  of  Q. 


=  1.  Points  shown  are  mapped  from  20 
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Fig.  19.  Response  of  elemental  filter  of  equation  (74). 
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Fig.  20.  Representation  of  frequency  sampling  filter  with  bandwidth  from 
300  cps  to  700  cps.  a)  comb  filter  b)  elemental  filter  c)  complete  filter. 
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Fig.  21.  Experimental  amplitude  vs  frequency  for  contiguous  band  pass 
filters  using  frequency  sampling  technique. 
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Fig.  22.  Digital  filter  with  quantization  error. 
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C62-569 


Fig.  24.  Canonical  representation  of  digital  filter  with  quantization  (drawn 
for  M  =  N). 
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